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Abstract 

We study the time evolution of a rotating, axisymmetric, viscous accretion flow around black holes 
using a grid based finite difference method. We use the Shakura-Sunyaev viscosity prescription. However, 
we compare with the results obtained when all the three independent components of the viscous stress 
are kept. We show that the centrifugal pressure supported shocks became weaker with the inclusion of 
viscosity. The shock is formed farther out when the viscosity is increased. When the viscosity is above a 
critical value, the shock disappears altogether and the fiow becomes subsonic and Keplerian everywhere 
except in a region close to the horizon, where it remains supersonic. We also find that as the viscosity 
is increased, the amount of outflowing matter in the wind is decreased to less than a percentage of the 
inflow matter. Since the post-shock region could act as a reservoir of hot electrons or the so-called 
'Compton cloud', the size of which changes with viscosity, the spectral properties are expected to depend 
on viscosity strongly: the harder states are dominated by low-angular momentum and the low viscosity 
flow with a significant outflows while the softer states are dominated by the high viscosity Keplerian flow 
having very little outflows. 

1 Introduction 

Pringle (1981) pointed out that in present of viscosity, most of the matter of the disc accretes into the black 
hole while most of the angular momentum is taken farther away by very Httle matter. Though this basic 
principle is the key to believe that the standard models of Shakura & Sunyaev (1973, hereafter SS73) and 
Novikov & Thorne (1973) are possible, there are very few numerical simulations which actually show how 
a Keplerian disc is formed, or whether it should always form. There are clearly two possibilities: when 
the companion supplies high angular matter (as in a low mass X-ray binary), viscous stress must be very 
large to transport this in order that the accretion may occur. On the other hand, when mainly winds from 
the companion (as in a high mass X-ray binary) is accreted the flow need not produce a Keplerian disc 
(Chakrabarti & Titarchuk, 1995, hereafter CT95; Smith, Heindl and Swank, 2002) since the viscosity may 
not be high enough to transport angular momentum rapidly to redistribute the angular momentum. It has 
been shown by Chakrabarti (1989) that low angular momentum, non dissipative flows produce axisymmetric 
standing shock waves in tens of Schwarzschild radii, but the presence of a large viscosity (Chakrabarti, 
1990, hereafter C90) will remove the shock wave since the Rankine-Hugoniot relation is not be satisfled in a 
highly dissipative flow. A conflrmation of such an assertion, originally made in the context of the isothermal 
flows, came both through numerical simulations (Chakrabarti & Moltcni, 1995, hereafter CM95) as well 
as theoretical studies of flows with more general equation of states (Chakrabarti, 1996a, hereafter C96; 
Chakrabarti & Das, 2004; Das & Chakrabarti, 2004). The theoretical works showed that the topology of the 
flow is modified when viscosity is added and beyond a critical viscosity, the centrifugal pressure supported 
shocks do not form. The numerical work of CM95, which was based on the smoothed particle hydrodynamics 
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(SPH) in a one dimensional flow, clearly demonstrated the transport of angular momentum by the viscosity. 
It showed that for low enough viscosity parameter a (SS73), the standing shock continues to form. However, 
with the increase in viscosity, the shock weakens and moves outward. Ultimately, for high enough viscosity 
the shock disappears. Most interestingly, the post-shock subsonic flow acquires a Keplerian distribution. 
Thus, for high enough viscosity, the entire disc becomes siibsonic and Keplerian, the assumed condition for 
the Shakura-Sunyaev disc. Igumenshchev et al. (1996) studied two-dimensional flows but concentrated only 
the inner region of the disk, namely, the region less than I^Vg. The main interest was to study the transonic 
nature of the flow just before the matter enters into the horizon. They find that a torus like structure forms 
close to the black hole and the angular momentum increases outward. Lanzafame et al. (1998) carried out 
the SPH simulation in a region with a radial extent of 50rg in two dimensions and concentrated on the shock 
formation. In particular, they showed that high viscosity can remove the shock waves from the accretion 
flows. However, SPH method is often intrinsically dissipative and it was necessary to verify the behaviour 
with other simulation methods. Igumenshchev et al. (1998) used the finite difference method and allowed 
the heat generated by viscosity heating to be radiated away or absorbed totally. The computational box was 
up to 300rg, but the outer boundary condition was that of a near Keplerian flow having no radial velocity. 
The inner boundary was kept at Sr^. Thus the possibility of having a shock or the inner sonic point were 
excluded. The disc was found to be stable for very high viscosity parameter a and less stable for lower 
a. Igumenshchev et al. (2000) extended earlier work by studying the dependence on the polytropic index 
7 which varied from 4/3 to 5/3 as well as viscosity parameter a and again found that the stability of the 
solutions depend on these parameters. 

So far, however, the fate of the solutions which include shocks have not been studied. The formation of 
shock required injection of very lower angular and low viscosity flows and the inner boundary must be inside 
2rg in order that the inner sonic point is formed (C89). In Giri et al. (2010, hereafter Paper I), the rcsidts 
of standing and oscillating shock formations in inviscid flows was presented. We now extend this study by 
adding viscosity. The theoretical work discussed in C89 or Chakrabarti (1996a) were carried out for a one 
dimensional flow since they required the sonic point analysis. For a two dimensional flow, a completely 
self-consistent theoretical solution is not possible. This is why a numerical simulation is necessary to answer 
the following questions: (a) Are the conclusions based on theoretical considerations continue to remain valid 
for a two dimensional flow? (b) Do the shocks survive for higher viscosity? (c) How does the outflow rate 
depend on viscosity? (d) Are all the components of the viscous stress important in a thick accretion flow? 
(e) Under what condition does the Keplerian flow form? 

In the next Section, we show how we introduce the viscosity into the grid-based flnite difference method 
(Molteni, Ryu & Chakrabarti, 1996). We consider both the cases, namely, when only the r<^ component is 
present, and when all the components, namely, r</), (^z and rz components are present. In §3 & 4, we present 
the methodology and the results of our simulation respectively. Finally, in §5, we draw concluding remarks. 



2 Model Equations with Viscosity 

In what follows, we consider a two dimensional axisymmetric flow around a Schwarzschild black hole. Instead 
of using general relativity, we use the well known pseudo-Newtonian potential prescribed by Paczynski 
and Wiita (1980). We use the cylindrical polar coordinates (r, (/> and z^. The mass, momentum and 
energy conservation equations in a compact form using non-dimensional units are given in Molteni, Ryu & 
Chakrabarti (1996; hereafter MRC96). We use the mass of the black hole Mbh (assumed to be ten solar 
mass in this paper), the velocity of light c and the Schwarzschild radius Vg = 2GMl(? as the units of the 
mass, velocity and distance respectively. 

The equations governing the inviscid flow have been presented in Ryu et al. (1995), MRC96 and Paper I 
in great detail and we do not repeat everything here. In the conservative form, the equations are given by, 
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Here, energy density E (without potential energy) is defined as, E = p/{'y—l)+p{v^ + Vg+v'^)/2, p is the 
mass density, 7 is the adiabatic index, p is the pressure, Vx, vq and Vz arc the radial, azimuthal and vertical 
component of velocity respectively. In the case of an axisymmetric flow without viscosity, the equation for 
azimuthal momentum states simply the conservation of specific angular momentum A, 

d\/dt = 0. 

In an inertial frame of reference, the general form of the equations of the flow (Batchelor 1967) is 
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where, v is the flow velocity, p is the fluid density, P is the pressure, r is the stress tensor, and Fb represents 
body forces (per unit volume) acting on the fluid and V is the Del operator. Typically body forces consist of 
only gravity forces, but may include other types (such as electromagnetic forces). Here t is the viscous stress 
having six mutually independent components. In cylindrical coordinates the components of the velocity 
vector given by v = {vr,v^,Vz)- The six independent components of the viscous stress tensor (Landau & 
Lifshitz 1959) are listed here in cylindrical coordinates, Trr,Tr,i,,TrzTT^^,T^z & Tzz- If we split all the viscous 
stress tensor, three components of equation 2 takes the following forms (Landua & Lifshitz 1959 and Acheson 
1990). The Vr component of Navier-Stokes equation is given by 
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Again, the component is given by 



dt ^ dr r r d(j) ^ dz 
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Finally, the Vz component is given by 
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where, /i is the dynamic viscosity defined hy ^ ~ rjp and rj is called the kinematic viscosity. Here, Fr,F^ 
& Fz are the so-called body forces. The only body force which is present in our system is the gravitational 
force. Thus, F^ = pG^, Fj- = pGr^F^ = pGz, where G^, Gr and Gz are the components of acceleration due 
to gravity, namely, 

1 r 

and 
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where, R = s/r"^ + z"^. For the axisymmetric case, = and thus, F^ = 0. As we have chosen the 
axisymmetric case, we have neglected ^ added terms. So the above equations reduces as following. 
Equation (3a) reduces to 

. dVr dVr dVr , 
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Equation (36) takes the form. 



Equation (3c) reduces to. 
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In case of a thin accretion flow, it is customary to use only Tr^ component since it is the dominant 
contributor to the viscous stress (SS73). This is responsible for transporting angular momentum along the 
radial direction. The other components are assumed negligible. Thus, considering only Tr^,, the only viscous 
term which goes into the Eq. (46) is 

dr'' r dr r^^ 

= ^■^(^''^r4>), (5) 
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where, 
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Here, f2 is angular velocity and defined as 
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It has long been suspected that the diffusion of angular momentum through an accretion flow is driven by 
turbulence. The 5'5'a-model (SS73) introduced a phenomenological shear stress into the equations of motion 
to model the effects of this turbulence. This shear stress is proportional to p, the total pressure. This shear 
stress permits an exchange of angular momentum between neighbouring layers. Using SSa prescription, we 
can assume r^^ = — ap, where a is a proportionality factor which need not be constant throughout the flow. 
This has been used in the simulations by CM95 and C96. In Eq. (5) we put —ap in the place of Trtj,. In that 
case, the viscous term reduces to 

4#(-V)=-a(^ + |^). (6) 

In case of thick accretion flow, all viscous stress could be signiflcant in flow dynamics. Assuming that 
the flow is thick /i r, the a-prescription could be written as (e.g., Igumenshchev et al. 1996) as 

M = oisPtt-, (7) 
where, as is constant of order 1, a is the adiabatic sound speed, and 

r or 

is the Keplerian angular velocity. We carried out simulations using this value of /x and either the rcf) 
component of the viscous stress or all the components of the viscous term. We clearly flnd major differences 
in our results. These will be discussed below. Note that, we used a different notation for the viscosity 
parameter to differentiate between cases where the Shakura-Sunyaev prescription (Eq. 6) is used or when 
the deflnition (7) is used for coefflcient of viscosity. 

It is to be noted that the heat generated by the viscous dissipation is assumed to be radiated away 
instantly. Thus we do not consider any effect of the heating on the dynamics of the flow in the present 
paper. 

In the literature, viscosity prescriptions other than those discussed here have also been tried out, especially 
when the shock is present. We followed also the prescription Macfadyen & Woosley (1999) where a was 
assumed to be constant when > Vr while it is assumed to be scaled as v^/vr in the pre-shock flow to 
reduce the shear. The result remains very similar to our present result. This is because the pre-shock flow, 
where > Vr, is cooler with a lower thermal prcssiirc and thus the angular momentum transport rate is 
weaker. In the post-shock flow, due to high thermal pressure, the angular momentum transport rate is higher 
and our disk becomes similar to Keplerian. Thus a constant a prescription plays a role similar to that in a 
Keplerian disk. 



3 Methodology of the Numerical Simulations 

The setup of our simulation is the same as presented in Paper I. We consider a viscoiis, axisymmetric flow in 
the Pseudo-Newtonian gravitational field of a point mass M(,h located at the centre in cylindrical coordinates 
[x, 9, z]. We assume that at infinity, the gas pressure is negligible and the energy per unit mass vanishes. We 
also assume that the gravitational field of the black hole can be described by Paczynski & Wiita (1980), 
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where, R = Vr^ + z^ and the Schwarzschild radius is given by 

Rg = 2GMbh/c'. 

We also assume a polytropic equation of state for the accreting (or, outflowing) matter, P = Kp'* , where, 
P and p are the isotropic pressure and the matter density respectively, 7 is the adiabatic index (assumed to 
be constant throughout the flow, and is related to the polytropic index n by j = 1 + 1/n) and K is related 
to the specific entropy of the flow s. In all our simulations, we take 7 = 4/3. 

The computational box occupies one quadrant of the r-z plane with < r < 50 and < z < 50. The 
incoming gas enters the box through the outer boundary (which is, in reality, of a vertical cylinder in shape), 
located at xi, = 50. We have chosen the density of the incoming gas pin = 1 for convenience since, in the 
absence of self gravity and cooling, the density is scaled out, rendering the simulation results valid for any 
accretion rate. We supply the radial velocity Vr, the sound speed a (i.e., temperature) of the flow at the 
boundary points. We take the boundary values of density from the standard vertical equilibrium solution 
(Chakrabarti, 1989). We inject matter at the outer boundary along the z-axis. We consider only the first 
quadrant and use the reflection boundary condition on the equatorial plane and z-axis to obtain the solution 
in other quadrants. In order to mimic the horizon of the black hole at one Schwarzschild radius, we placed an 
absorbing inner boundary at R = l-lrg, inside which all material is completely absorbed into the black hole. 
For the background material, we consider of density pbg = 10~^ having the sound speed (or, temperature) 
to be same as that of the incoming matter. Hence, the incoming matter has a pressure 10^ times larger than 
that of the background matter. This matter is put to avoid unphysical singularities caused by 'division by 
zero'. This is washed out and replaced by the incoming matter withing a single dynamical time scale. The 
calculations were performed with 512 x 512 cells, so each grid has a size of 0.097 in units of the Schwarzschild 
radius. The timescale of matter from the outer boundary is about 0.07s as computed from the sum of 
dr/ <Vr> over the entire radial grid, <Vr> being averaged over 20 vertical grids. 

All the simulations have been carried out assuming a stellar mass black hole (M = IOMq). We carry out 
the simulations for several hundreds of dynamical time-scales. In reality, our simulation time corresponds to 
a few seconds in physical units. The conversion of out time unit to physical unit is 2GM/(?, and thus the 
physical time for which the programme was run would scale with the mass of the black hole. 

4 Simulation Results 

The total variation diminishing (TVD) method that we use here was primarily developed to deal with fluid 
dynamics (Harten, 1983). In the astrophysical context, Ryu et al. (1994 & 1996) developed the TVD scheme 
to study astrophysical flows around black holes. In Paper I, the oscillation phenomena in accretion flows 
around black holes related to the QPO were reported. That study of an inviscid accretion flow around 
black hole showed that the shock location changes with the change of specific angular momentum (A) and 
specific energy {£), both of which were constant. In the present situation of a viscous flow, none of these is 
constant. Recently, one dimensional solution for quasi-spherical viscous flow was investigated by Lee, Ryu 
& Chattopadhyay (2011) who found that there are oscillatory propagating shocks moving outward when the 
viscosity is large enough. 

4.1 Isothermal injection at the outer boundary 

As mentioned earlier, we chose the outer boundary of the simulation grid at r = 50. The speciflc angular 
momentum (A) of the flow is chosen to be 1.66 (For comparison, we note that the marginally stable angular 
momentum is 1.83 in this unit.) and the specific energy {£) of the flow at the equatorial plane (z = 0) 
is chosen to be 0.0035. The radial velocity pointing to the origin is chosen to be constant in all heights 
V = (vr^ + Vz^)^ = 0.072. These injection parameters correspond to those of low angular momentum 
transonic flow solutions and are qualitatively different from those of earlier workers such as Igumenshchev et 
al. (1998, 2000) since they mainly inject near Keplerian flow with no radial velocity. From the energy of the 
flow, we obtain the sound speed at the equatorial plane to be 0.059. We employ an isothermal outer boundary 
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Figure 1: Variation of (a) Mach nuniber and (b) density with radial distance in a viscous transonic flow on 
the equatorial plane. As the viscosity parameter is increased, the angular momentum is transported outward 
shifting along with it the centrifugal pressure supported shock wave. The a parameters are [(left to right in 
(a) and bottom to top in (b)]: 0.0, 0.0175, 0.035, 0.0525, 0.06125, 0.07, 0.0735 respectively. 

condition (MRC96). In other words, we take the same sound speed at all heights at the outer boundary. We 
add the SSa viscous term in the non- viscous system as discussed earlier. Wc stop the simulation at t = 24.75 
seconds. This is more than three hundred times the dynamical time. Thus, the solution has most certainly 
come out of the transient regime and started exhibiting solutions characteristics of its flow parameters. The 
simulation results will be discussed now. 

In Figs, la and lb, we compare the Mach number and the specific angular momentum variations in the 
equatorial plane of the flow for various a. To make the comparison meaningful, all the runs were carried out 
up to t — 24.75s. Each result is obtained starting with an inviscid flow a — (marked) and then gradually 
increasing a till the shock goes out of the grid and eventually disappears. The values of a for which the 
curves are drawn are (left to right in Fig. la and bottom to top in Fig. lb): 0.0, 0.0175, 0.035, 0.0525, 
0.06125, 0.07, 0.0735 respectively. The shock location shifts outward with viscosity exactly as predicted in 
COO and C96. For an inviscid flow (a = 0), the matter bounces back from the centrifugal barrier and flows 
outward near the equatorial plane. As viscosity is enhanced, the angular momentum is transported outward 
and hence specific angular goes up. The velocity goes down and density goes up. This can be seen in Fig. 
lb in the 3 — 30 region. In this particular example, the shock disappears above ~ a = 0.074, which is the 
critical value of a here. Since the grid boundary is at a finite distance in a numerical grid it is difficult of 
show the disappearance of the shock since that would require changing the boundary condition dynamically 
to let the shock through when it reached the boundary. The results are indeed similar to the results of CM95 
where one dimensional viscous isothermal fiow was treated using Smoothed Particle Hydrodynamics (SPH) 
code. As described in C90, the critical a (ac) is defined by that specific a. for which the subsonic branch of 
the transonic fiow solution passes through both the inner and outer sonic points. Thus ac clearly depends 
on the injected flow parameters. Hence, for a diff'erent choice of injected parameters, the critical value of a 
will be different. 

In Fig. 2, we show how the density of matter and the velocities vary with viscosity. We superpose 
contours of constant density and velocity. The length of the arrows are proportional to velocity, the longest 
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being that oi v = 0.6. The results arc for t = 24.75s. a = 0.0 (a: top loft), 0.0525 (b: top right), 0.0735 
(c: bottom left) 0.074 (d: bottom right) respectively. The contour having minimum density has pmin = 0.2. 
The maximum density is (a) pmax = 450, (b) Pmax = 395, (c) pmax = 354 and (d) Pmax = 286. The contour 
interval is 6p = 0.25. We observe that the standing shock is truly two dimensional and oblique (Chakrabarti, 
1996b). Here three shocks meet at a point, is a good example of a prominent triple-shock which forms away 
from the equatorial plane. However, it is not a pure triple-shock, since the matter flows from right to the 
left for both the shocks facing the upstream. As the matter flows in from the right hand side and passes 
through the shock, it becomes hot and puffs up as a thick accretion disc (see also, Molteni, Lanzafame and 
Chakrabarti, 1994, hereafter MLC94). The standing shock moves outwards as the viscosity parameter is 
increased. For a = 0.074 which is larger than the critical a for this case and the shock disappears. In a 
numerical simulation with a finite injection radius the shock, albeit very weak, stays very close to the outer 
boundary. The piled up angular momentum in the post-shock fiow clearly drives the shock outwards. In 
Fig. 3 we plot the time variation of the shock location at the equatorial plane for various viscous parameter 
a. The values of a are (from bottom to top): 0.0, 0.0175, 0.035, 0.0525, 0.06125, 0.0735 respectively. It is 
clear that the average shock locations are shifting outwards when the value of a increase. Shocks close to the 
black hole exhibit lower amplitude and higher frequency oscillations while those farther out show opposite 
effects. This is mainly because the frequency is decided by the inverse of the irifall time in the post-shock 
region. The compression wave in the post-shock region bounces back from the centrifugal barrier and pushes 
the shock outward. At some point the outward journey is stopped when post-shock pressure drops and the 
shock turns back. Most interestingly, for a = 0.0525 and a = 0.06125, the oscillation of the shock disappears 
and standing shocks are formed while for other as there are oscillations. This is not surprising, since, as 
was shown for the inviscid ffow (Chakrabarti, 1989) as well as the viscous flow Chakrabarti & Das (2004), 
the Rankine-Hugoniot relations are satisficid only in a limited region of the parameter space, and beyond 
the critical viscosity the relation not satisfied at all. In this case, the flow has two sonic points and the 
high entropy solution must pass through the inner sonic point (C89). Hence, the fiow generates entropy 
through a shock jump, but the location itself cannot be fixed because the Rankine-Hugoniot condition is not 
satisfied. This creates an unstable situation and consequently, the shocks can oscillate (see, Ryu, Molteni & 
Chakrabarti, 1996; Paper I). 

It is interesting to study the effects of viscosity on the outflow rate. We observed that with viscosity 
the shock recedes and becomes weaker. In Chakrabarti (1999), it was suggested that the ratio Rm of the 
outflow to the inflow rate would be guided by the compression ratio at the shock. In Fig. 4, we plot the 
time variation of the ratio of the outflow to inflow rate as the viscosity is enhanced. The values are same as 
in Fig. 3 (from the top curve to the bottom curve). We clearly notice that although the ratios exhibit short 
time scale fluctuations, the average values decrease as viscosity is enhanced. This could have been guessed 
also from Fig. 2(a-d), where the lengths of the outgoing arrows in Figs. 2(a-b) are reduced in number and 
size in Figs. 2(c-d). The longest arrow corresponds to v ^ 0.6. 

4.2 Vertical equilibrium at the outer boundeiry 

We now use a different model for injection in order to impress that the basic results remain the same even 

when we assume the flow to be in vertical equilibrium (Chakrabarti, 1989) at the outer boundary. We inject 
through 50 grids out of 512 grids, i.e., the height of the disc at injection is nearly 5 Schwarzschild radii. In 
this case, the injection rate of the momentum density is kept uniform throughout the injected height at the 
outer edge. The specific angular momentum (A) and specific energy (E) at the outer boundary is chosen to 
be the same as in the previous case. Here too we stop the simulation a,t t = 24.75 second. The results of the 
simulation are discussed now. 

In Figs. 5a and 5b, we show the distribution of Mach number on the equatorial plane and that averaged 
over 15 grids from the equatorial plane when a is 0.0, 0.018, 0.0225, 0.0315, 0.05 and 0.09 respectively 
(from leftmost curve to the rightmost curve). As the a-parameter is increased, the fiow behaviour changes 
dramatically. As in the the previous case, the shock rapidly propagates outward due to the faster transport 
of angular momentum in the post-shock flow as compared to the pre-shock flow. However, since in the 
vertical equilibrium model, the ram pressure of the injected flow is less, the turbulence and the back flow on 
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Figure 2: Changes in the density and velocity distribution with the change of viscous parameter a at 
t = 24.75s. Here, densities are in normalized unit, radius and velocity are in Schwarzschild unit. Here, 
a = 0.0 (top left), 0.0525 (top right), 0.0735 (bottom left) 0.074 (bottom right)respectively. For details see 
the text. 
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Figure 3: Variation of the shock location at the equatorial plane with time for various viscous parameter a 
(from bottom to top: 0.0, 0.0175, 0.035, 0.0525, 0.06125, 0.0735 respectively). 
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Figure 4: Time variation of the ratio of the outflow to inflow rates as viscosity parameter is enhanced. 
Though there are short time scale fluctuations, the average values decrease as viscosity is increased showing 
a direct relation of the outflow rate with the strength of the shock. The viscosities are same as in Fig. 3 
(from the top curve to the bottom curve). 
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Figure 5: Variation of the (a) Mach number on the equatorial plane and (b) the Mach number averaged 
over 15 vertical grids as a function of the radial distance on the equatorial plane as the a is increased from 
(leftmost curve) to 0.09 (rightmost curve). See the text for details. 

the equatorial plane remains important even for high a. 

Figures 6(a-d) show the distribution of average specific angular momentum of the flow (density weighted 
average over 15 vertical grids from the equatorial plane) at the end of our simulation for different values of 
a (thin solid curve). For comparison, we plot the specific angular momentum distributions of a Keplerian 
orbit (thick upper curve). We also plot a 'Keplerian' distribution (thick lower curve) at a height of 200 grids 
(^ 20 Schwarzchild radii) and compare the 15 vertical grid average angular momentum distribution that is 
obtained from the simulation at that height (dashed curve). The latter 'Keplerian' distribution was obtained 
by equating the horizontal component of the gravitational force at that height with the centrifugal force. 
The results are shown at i = 24.75s. The a-paramcters chosen for the Figs. 6a, 6b, 6c and 6d arc 0.0, 0.01, 
0.05, and 0.09 respectively. We note that as the viscosity is increased, the distribution in the post-shock 
region gradually becomes closer to the Keplerian distribution, although, below r = 3, the distribution is 
always highly sub-Keplerian. At the shock, the distribution shows a jump. This is because for a given a, the 
transport rates in the pre- and the post-shock flows are different, being very high in the post-shock region 
due to higher pressure. For the high enough viscosity, when the shock reaches infinity (large distance), the 
angular momentum distribution becomes that of a Keplerian flow, which is what is expected. This is thus 
one scenario by which a Keplerian disc may form in a highly viscous flow. 

4.3 Viscous flow with a constant height injection in radial direction 

For the sake of comparison with the results of inviscid flows presented in Paper I, we simulate the case 
where the injection is purely in the ~X direction. The injected flow at the outer boundary has the same 
sound speed (temperature) as obtained from the theoretical 'constant height' model. The specific energy 
and specific angular momentum arc E ~ 0.003 and A = 1.76 respectively. The simulation was carried out 
up to t = 7.63s, or about one hundred dynamical time (computed as a sum of dr/ < Vr > over the whole 
radial grid, < Vr > being averaged over 20 vertical grids). 

Figure 7 shows the distribution of Mach number along the equatorial plane for a = 0.0 (dashed), 0.02 
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Figure 6: Comparison of the specific angular momentum distributions (thin curves) of an accretion flow as 
viscosity parameter a is varied, a = (a) 0.0, (b) 0.01, (c) 0.05 and (d) 0.09. For all the cases, the result is 
compared with the Keplerian angular momentum distribution (thick curves). 
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Figure 7: Radial distribution of the Mach number on the equatorial plane for a flow with a constant height 
injection in the radial direction at the outer boundary. The viscosity parameters are: a = 0.0 (dashed), 0.02 
(solid), 0.035 (long dashed), 0.07 (dash-dotted) and 0.09 (log dash-dotted). The inner shock disappears at 
about a. ~ 0.05 and the outer shock disappears at a ~ 0.1. 



(solid), 0.035 (long dashed), 0.07 (dash-dotted) and 0.09 (log dash-dotted) respectively. As in Paper I, we 

see the formation of a strong (outer) shock at ^ 25 and a weak (inner) shock at ^ 3 when the flow is inviscid. 
For a = 0.02 and 0.035 both the shocks still form, albeit being farther out and weaker. For higher viscosity, 
the inner shock completely disappears. The outer shock disappears at even higher viscosity. This shows that 
the critical viscosity for the inner shock is much lower than that for the outer shock. Since, theoretically 
only the outer shock was predicted, this is therefore a totally new result and could not have been anticipated 
without numerical simulations. In Fig. 8, we show the density and velocity distributions for a = (top left) 
and 0.02 (top right) respectively at the end of the simulation. We zoom the inner region to show how the 
inner shock has shifted from ~ 3 (bottom left) to ~ 8 (bottom right) as viscosity is increased. 

4.4 Effects of boundary location 

It may be noted that we have run the simulations above for viscosity parameters slightly below the critical 
viscosity, since the shock is nearing the outer boundary. As we have injection of matter at the outer boundary, 
it would prevent the shock from leaving the grid near the injection area. In order to prove that the shock 
actually runs away and disappears, one requires an infinitely large bomidary. However, that this must happen 
can be easily shown by running a few cases using the viscosity parameter from both sides of the critical value. 
In Fig. 9(a-b), we present two such cases with (a) E = 0.0035, A = 1.66, Xh = 50 and (b) £ = 0.002, A = 1.66, 
Xb — 100. The viscosity parameters are marked on the curves. The critical viscosities are etc ~ 0.0738 and 
0.0325 respectively. We clearly see that for a < Uc, the shock first goes out farther before returning back 
and settling down at a certain finite distance with some small amplitude oscillation. However, for a > Uc, 



14 



Figure 8: The density and velocity distributions for a = (top left) and 0.02 (top right) respectively at the 
end of the simulation where the flow was injected with a constant height. We zoom the inner region to show 
the shifting of the inner shock from 3 (bottom left) to 8 (bottom right) when viscosity is increased. 
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the shock never returns back and continues to go outward. This behaviour is independent of the location of 
the boundary. 

4.5 Effects of viscous stress components 

In this case, we used Eq. (7) in Eqs. 4(a-c). We considered these cases to understand the importance of 
various components of the viscous stress. In Fig. lO(a-d), radial distribution of the Mach number and the 
specific angular momentum on the equatorial plane and that averaged over fifteen vertical grids from the 
equatorial plane are shown. The dashed-dotted curve shows the results when only the rcj) component is used 
and the solid curve shows the results when all the three components are used. The results are plotted at 
t = 10s. The fiow parameters arc £ = 0.035, A = 1.66 and ag ~ 0.03. We also plot the results for the 
inviscid flows for comparison (dotted curve). Note that the distribution of angular momentum inside the 
shock remains almost constant when a single (^np) component is used while it becomes similar to Keplerian 
when all the three components are included. 

In Figs, ll(a-d), we plot the distribution of density and velocity when only rcj) component of the viscous 
stress is used (top left). Note that the jaggedness of the shock goes away when all three components of 
the viscous stress are included (top right). A density maximum occurs in the post-shock region. Thus the 
post-shock region behaves like a thick accretion disk (Paczyhski & Wiita, 1980) as was also pointed out in 
SPH simulations (MLC94). The post-shock region, formed purely due to the centrifugal force, is known as 
the CENtrifugal pressure dominated BOundary Layer or CENBOL. This region is believed to be responsible 
to emit hard X-rays in black hole candidates by inverse Comptonizing soft photons coming from Keplerian 
disks believed to form near the equatorial plane where the viscosity is high (CT95). The corresponding 
specific angular momentum distributions are in the bottom-left and the bottom-right panels respectively. 
Note that up to the disk center the angular momentum rises rapidly and becomes almost Keplerian as shown 
in Fig. ll(c-d). The parameters chosen are the same as in Fig. 10. A comparison of the behaviour with 
and without all the three components indicates that the behaviour becomes smoother with transport of 
momentum takes place in all directions. The outflow also has a larger specific angular momentum. 

5 Discussions 

In this paper, we have presented the results of the numerical simulations of a two dimensional, axisymmetric, 
viscous accretion flows. Three parameters, namely, the specific angular momentum, the specific energy and 
the viscosity parameter dctcirminc the complete solution, although results depend somewhat on the injection 
process at the outer boundary. While both inviscid (Chakrabarti, 1989) and viscous (Chakrabarti & Das, 
2004) flows allow solutions with or without centrifugal barrier dominated shock waves, we concentrated 
mostly on the cases when the shocks arc formed. Wc find that the shocks move outward as the viscosity is 
enhanced and the post-shock region roughly attains a Keplerian distribution. When the viscosity parameter 
is very high, the shock moves to a large distance and the whole disc becomes a Keplerian disc. We also 
found that the condition of standing shock wave formation may be satisfied only in a narrow range of the 
viscous parameter (keeping other parameters as constants) which is in line with the conclusions drawn in 
Chakrabarti & Das (2004). When the Rankine-Hugoniot conditions are not satisfied, the shocks tend to 
oscillate (Ryu, Chakrabarti, Molteni, 1997; Paper I) and the frequency of oscillation is decreased and the 
amplitude is increased as the shock moves out. 

In contrast to these simulations, the simulations of Igumenshchev and collaborators (1996, 1998, 2000) 
concentrated on shock-free solutions. Either the simulation boundary was too small or too big with high 
viscosity, the shocks did not form in these simulations. Thus a direct comparison is not possible. In some 
of the models (e.g., CT95) the post-shock region or CENBOL is treated as the Compton cloud which plays 
a major role in shaping the spectrum of the accretion fiows. Similarly, the CENBOL is thought to be 
responsible for the outflows. In Chakrabarti (1999) it was estimated that the ratio of the outflow to the 
inflow varies with the strength of the shock. In the present case, we varied the strength of the shock by 
increasing viscosity and found that the outflow rate, though fluctuating in small time scales, becomes smaller 
as the viscosity parameter is increased. 
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Figure 9: Evolution of the shock location as the viscosity parameter a is changed for two different boundary 
conditions: (a) xt = 50,£ = 0.0035, A = 1.66; (b) Xb = 100,5 = 0.002, A = 1.66. The critical viscosities are 
Qc = 0.0738 and 0.0325 respectively in these cases. In both the cases, when a. < a,., the shock first goes 
out farther before returning back and settling down, while for a > ac, the shock never returns back and 
continues to go outward. 
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Figure 10: Radial distribution of the (a & b) Mach number, and the (c & d) specific angular momentum on 
the equatorial plane (a & c) and that averaged over fifteen vertical grids from the equatorial plane (b& d) 
for the cases when only the rcj) component (dashed-dotted) and all the three components (solid) are used. 
The results arc plotted at t = 10s. We also plot the results for the inviscid flow for comparison (dotted). 
The angular momentum inside the shock remains almost constant when a single (r0) component is present 
while it becomes similar to Keplerian when all the three components are included. 
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Figure 11: Distribution of density and velocity when only rt^ component of the viscous stress is used (top 
left). Note that the jaggcdncss of the shock goes away when all three components are included (top right). 
Density maximum and a consequent thick accretion disk like structure is formed in the post-shock region in 
the latter case. The corresponding specific angular momentum distributions are in bottom left and bottom 
right respectively. The parameters are same as in Fig. 10. 
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In Paper I, we emphasized the formation of a weaker inner shock eloscr to the black hole. This was not 
predicted by the theoretical works. In the present paper, we find that this inner shock also becomes weaker 
and moves outward when the viscosity is introduced. The critical viscosity parameter for removal of the 
inner shock is lower than what is needed to remove the outer shock. 

In Molteni, Sponholz and Chakrabarti (1996) and Chakrabarti, Acharyya & Molteni (2004) it was shown 
that the shocks oscillated when the cooling time scale matches with the infall time scale. The simulation 
did not include viscosity. In Ryu, Chakrabarti and Molteni (1997) no cooling or viscosity was added, but 
the shocks were shown to oscillate when the Rankine-Hugoniot relation was not satisfied. In the present 
paper, for the first time, we show that even for a viscous flow, the same conclusion holds true as long as the 
viscosity is less than the critical value and the shocks are not totally removed. 

It is well known that the light curves of the radiation coming from an accretion disc around a black hole 
exhibit quasi-periodic oscillations (QPOs). It appears that in the absence of a hard surface, the standing 
shock itself behaves like a hard surface and its oscillation changes the size of the post-shock region signifi- 
cantly. This could be the cause of the low and intermediate frequency quasi-periodic oscillations (Chakrabarti 
& Manickam, 2000; Rodriguez et al. 2004; Remillard et al. 2006; Gliozzi et al. 2010; Qu et al. 2010). When 
the viscosity is increased, the Keplerian rate is enhanced and at the same time, the shock recedes to a large 
distance and the time period is increased. Asymptotically, this means that a Keplerian disc should not show 
QPOs. We also observe that away from the equatorial plane the angular momentum is sub-Keplerian (i.e., 
smaller compared to the specific angular momentum on the equatorial plane). A picture similar to this 
has been postulated in the two component advective flow model (GT95) where sub-Keplerian flows flank a 
Keplerian disc on a equatorial plane. However, a complete picture would emerge only after radiative transfer 
is added. 

The authors thank D. Ryu and I. Chattopadhyay for providing helps to include viscosity in the original 
code for non-viscous flow used in Paper I. They are also thankful to an anonymous referee for helpful 
suggestions. 
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